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Introduction 

The importance of reducing jet noise in both commercial and military aircraft applications has 
made jet acoustics a significant area of research [1]. A technique for jet noise prediction commonly 
employed in practice is the MGB approach [2], based on the Lighthill acoustic analogy [3], This 
technique requires as aerodynamic input mean flow quantities and turbulence quantities like the ki- 
netic energy and the dissipation. The purpose of the present paper is to assess existing capabilities for 
predicting these aerodynamic inputs. Two modem Navier-Stokes flow solvers, coupled with several 
modem turbulence models, are evaluated by comparison with experiment for their ability to predict 
mean flow properties in a supersonic jet plume. Potential weaknesses are identified for further inves- 
tigation. Another comparison with similar intent is discussed by Barber et al. [4]. The ultimate goal 
of this research is to develop a reliable flow solver applicable to the low-noise, propulsion-efficient, 
nozzle exhaust systems being developed in NASA focused programs. These programs address a 
broad range of complex nozzle geometries operating in high temperature, compressible, flows. 

Seiner et al. [5] previously discussed the jet configuration examined here. This convergent- 
divergent nozzle with an exit diameter of 3.6 in. was designed for an exhaust Mach number of 
2.0 and a total temperature of 1680° F. The acoustic and aerodynamic data reported by Seiner et 
al. [5] covered a range of jet total temperatures from 104° F to 2200° F at the fully-expanded 
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nozzle pressure ratio. The aerodynamic data included centerline mean velocity and total temperature 
profiles. 

Computations were performed independently with two CFD codes, ISAAC [6] and PAB3D [7]. 
Turbulence models employed include the k-e model [8], the Gatski-Speziale algebraic-stress model 
[9] and the Girimaji model [10], with and without the Sarkar compressibility correction [11]. Cen- 
terline values of mean velocity and mean temperature are compared with experimental data. 

Codes and Models 

ISAAC (Integrated Solution Algorithm for Arbitrary Configurations) [6] is a finite-volume code 
of second-order accuracy which solves the full Favre-averaged Navier-Stokes equations. An upwind 
scheme based on Roe’s flux-splitting is used for the convective terms, central differencing for the 
diffusion terms and an implicit, spatially split, approximate-factorization scheme for iteration. 

The PAB3D code [7] solves the Reynolds-averaged Navier-Stokes equations in a thin-shear-layer 
approximation. Diffusion terms are central-differenced. The implicit iteration operator employs the 
van Leer scheme and the explicit terms ( e.g ., convective terms) are evaluated with the Roe scheme. 

The grid used here is composed of five blocks, one in the nozzle (61 x 61 mesh points), one 
external to the nozzle (61 x 61 mesh points) and three blocks downstream of the nozzle exit (65 x 121, 
97 x 121 and 97 x 121 mesh points). The mesh is very fine near the nozzle walls and near the jet 
axis in the plume, and becomes gradually coarser away from the axis. Grid-independence tests 
(see below) indicated this grid was adequate for the computations of the quantities examined in 
this work. Velocity, pressure and temperature boundary conditions consistent with experimental 
conditions are imposed at the upstream boundary inside the nozzle and a modest free-stream flow of 
M = 0.05 is imposed through a one-dimensional characteristic far-field boundary condition at all 
non-wall boundaries except the downstream boundary, where a subsonic outflow boundary condition 
is imposed. No-slip boundary conditions are imposed at all walls and axisymmetry is assumed. 

The turbulence models employed in this investigation include both k-e and algebraic-stress 
two-equation models. The Reynolds stress, r^, is modeled by the expression 

pTij = \pk5ij - ypCpfn — [{Sa - -Skk^ij) + <*4 —{SikWkj + SjkWki) 
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Sij = (l/2)(uij + Uj t i) andWij = (l/2)(u lJ — Uj ti ) being the symmetric and anti-symmetric parts of 
the mean-velocity-gradient tensor. Commas denote differentiation, the w* are velocity components 
and p is the average density. The remaining symbols are defined below. 

The turbulent kinetic energy k and the solenoidal part of the dissipation rate e s (the total dis- 
sipation e = e s + e c is the sum of the solenoidal and curl-free dissipations; in the case of no 
compressibility correction £ c = 0) are determined using the conservation equations 

~b {P^j^),j = P^mn^m,n P& "b ( (p pC^ffi J 
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For the k - £ model, a 4 = 0, a 5 = 0 and C* = = const. In the k and £ s equations, a k , a £ , 

C £l and C e 2 are constants and integration to the wall is possible using the damping functions /^, /i 
and / 2 . Also in the £ equation, we have i — £ s — p/p \V\/k\ 2 and \w , a wall-correction function. 

The Gatski-Speziale algebraic-stress model [9] employs (1) with C*, a , and a 5 functions of k, 
£ and invariants of Sij and Wij. No wall damping of the Reynolds stress is required, so / M = 0, and, 
in ISAAC, new f x and / 2 are introduced to preserve proper log-law behavior. 

The Girimaji algebraic-stress model [10] also employs (1); in this case the coefficients C*, a 4 
and tt 5 are determined via the (explicit) solution of a cubic algebraic equation; see the work of 
Girimaji [10] for details of this solution and of the selection of the physically relevant root. 

The models were tested both with and without the Sarkar compressibility correction [11], which 
gives for the compressible part of the turbulent dissipation £ c the value a c M?£ s , M t being the tur- 
bulent Mach number s /ri'i/a and a being a constant taken here to be 0.5. 

Results and Discussion 

Computations were performed at jet total temperatures of 104° F and 1550° F. Extensive tests of 
the adequacy of the spatial resolution were conducted; comparison of runs with the full grid and with 
the grid coarsened by factors of 2 and 3 showed negligible differences. The value of y + for points 
adjacent to the walls of the nozzle was always less than 0.5 for the full-resolution computation, 
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indicating that the wall layers were satisfactorily resolved. All runs were continued well past the 
point where a steady state in the mean quantities had been reached. 

Data for the 1550° F case appear in Figure 1 . The mean velocity and temperature at the centerline 
are shown as computed using the various models discussed in the previous section and compared 
with experimental results [5]. While there is a good deal of variation in the computational results 
(this variation is much greater than the experimental uncertainty [5]), the ISAAC results are better, 
overall, than the PAB3D results. The centerline mean velocity for the 104° F case is shown (model 
by model) in Figures 2 and 3 and agreement is fairly good for both codes. While the results shown 
here are worse than one would like, they are better than earlier results [4]. 

The simulations employing the Sarkar compressibility correction were generally no better than 
those not. In the 104° F jet case, the correction makes the PAB3D predictions significantly worse, 
and the ISAAC predictions are slightly high by about as much as they were low for the computa- 
tions without the compressibility correction. In the 1550° F jet case, the compressibility correction 
significantly improves the PAB3D results, but makes the ISAAC results dramatically worse. 

Figure 1 shows a grouping of the results from each of the two codes; this grouping occurred 
in the other case as well. It suggests the differences in results may be due more to differences in 
the codes than to differences in the models. In order to further investigate these differences, the 
k — e model in PAB3D was implemented in ISAAC (except for the PAB3D form for which 
involves derivatives not readily available in the ISAAC code), as was the PAB3D implementation 
of the Gatski-Speziale ASM. Comparisons between the ISAAC implementations of the models, the 
PAB3D implementations of the models and the PAB3D models implemented in ISAAC are shown 
in Figures 2 and 3 for the 104° F case. The results from ISAAC with its own models and with 
the PAB3D models were very similar for the k - e model and the Gatski-Speziale algebraic-stress 
model, indicating that the differences between the ISAAC and PAB3D results are primarily due to 
the codes themselves. 

In summary, overall agreement of the computations with experiment is good. The two codes 
each gave fairly consistent results with the different turbulence models, and the differences between 
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the codes seemed to be greater than the differences between the models. Additional evidence for this 
was given by the computations with ISAAC using the PAB3D versions of the models, which gave 
results much closer to the ISAAC results with its own models than to the PAB3D results. Possible 
reasons for these differences between the results of the two codes include different handling of 
viscous fluxes (thin shear layer in PAB3D vs. full Navier-Stokes in ISAAC), first-order advection 
in the PAB3D turbulence equations vs. second-order advection in all equations in ISAAC and other 
differences in the numerical algorithms employed in the two codes. 

Finally, we note that the aerodynamic input is only part of the story, and the MGB noise predic- 
tion may emphasize or de-emphasize different aspects of the input error. 
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K. S. Abdol-Hamid with PAB3D is gratefully acknowledged. P. Pao and J. Carlson supplied the 
grid. The ISAAC and PAB3D CFD codes were made available by NASA Langley Research Center. 
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ABSTRACT 


Through the use of both an a priori analysis of direct numerical simulation data and 
experiments with large-eddy simulations, a non-Smagorinsky grid-size dependence is estab- 
lished for the Smagorinsky sub-grid scale model for low numbers of resolved scales. It is 
shown that an increase in the Smagorinsky constant as the grid size is increased permits suc- 
cessful large-eddy simulations for cut offs approaching the energy-containing range of length 
scales. A detailed comparative analysis is made of the second-order turbulence quantities as 
determined by the differently resolutioned large-eddy simulations. 
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is periodic and shares many properties with common shear flows, but lacks many of the 
computational difficulties presented by homogeneous shear flow or such wall-bounded shear 
flows as channel flow or boundary-layer flow. Kolmogorov flow was originally proposed by 
Kolmogorov to his students as a simple flow for the study of stability problems in shear flows 
and has since been subjected to a number of stability analyses [5]. Turbulent Kolmogorov 
flow in two dimensions was studied by She [6]; three-dimensional direct numerical simulations 
of turbulent Kolmogorov flow employing the physical viscous stresses were performed by 
Shebalin and Woodruff [7]. Three-dimensional hyperviscosity simulations were performed 
by Borue and Orszag [8]. Some preliminary results from the present investigation have been 
presented earlier [9] and the present results have been recently presented [10]. 

The question of possible alternative grid-size dependences for the Smagorinsky model is 
here approached in two ways. First, an a priori analysis of a 64 3 direct-numerical simulation 
(DNS) data set is conducted. The fully-resolved velocity field from the DNS is filtered over 
a wide range of filter widths; for each filter width, both the Smagorinsky model formula 
and the Reynolds stress it is supposed to predict are evaluated. Comparing the two permits 
the filter width (or grid size) dependence of the Smagorinsky model to be assessed. The 
second approach to assessing the grid-size dependence of the Smagorinsky model is by LES 
experimentation: at each of a number of numerical resolutions, simulations are run with 
different values of the Smagorinsky constant. Determining which value of the Smagorinsky 
constant leads to the best results at a given resolution yields an empirically determined 
grid-size dependence. 

Encouragingly, both approaches lead to similar conclusions about grid-size dependence. 
In both cases, an intermediate range of wave numbers (or grid sizes) was found where the 
standard A 2 grid-size dependence worked well. Beyond the low-wavenumber (coarse res- 
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olution) end of this range, an enhanced grid-size dependence is found to be necessary to 
give the best possible results, both in the a priori analysis of the DNS data set and for 
the LES experiments. Satisfactory LES results were achieved at suprisingly low resolutions 
when an appropriate grid-size dependence was used, though there was deterioration in the 
predictions, particularly of the turbulent shear stress. 

The investigation outlined in this paper has sought to establish in a simple way that 
there is merit to the idea of altering the grid-size dependence of a sub-grid-scale model 
(in this case, the Smagorinsky model, with its traditional grid-size dependence based on 
inertial-range dynamics) in order to compensate for limited under-resolution in a large-eddy 
simulation. In keeping with this goal, a simple criterion for selecting the best Smagorinsky 
constant has been chosen (comparing the kinetic energies in different simulations) and sim- 
ple criteria for assessing the adequacy of the under-resolved simulations (comparing several 
second-order statistical quantities) were employed in this preliminary investigation. More 
sophisticated approaches to fixing the Smagorinsky constant (including a theoretical discus- 
sion of its grid-size dependence as the energy-containing range is approached from above) 
and to assessing the usefulness of under-resolved large-eddy simulations are the subject of 
current investigation. 

Following the problem statement and notational definitions of the next section, the a 
priori analysis of the DNS data set and the LES experiments are discussed. The paper 
concludes with a discussion of the results and their implication for the LES of complex 
flows. 

Problem Statement and Numerical Considerations. 

Kolmogorov flow is described by the Navier-Stokes equations subject to an artificial 
periodic body force with wavenumber k f and amplitude kjv 2 0 . An appropriate nondimen- 
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sionalization consists of using 1 fkj for the reference length and \/k } v 0 for the reference 
time; the Reynolds number is then v 0 /kjv. Taking the force vector to be in the positive x 
direction, the Navier-Stokes equations for an incompressible flow become 

Dv 1_ 1 _o ? . ,,\ 

— = — Vp + — V 2 v + i snu, (1) 

Dt p Re 

Periodic boundary conditions are imposed in all three directions, so the forcing wavenumber 
kf reappears in the formulation in the ratio of the forcing wavenumber to the wavenumber 
of the longest wave contained in the box; this ratio will hereafter be denoted simply k s . 

The simulations examined in the present investigation have a Reynolds number of 28 and 
a forcing wavenumber of 6. This Reynolds number is sufficiently high to give fully developed 
turbulence (transition to turbulence occurs at Re = 12 - 13 in this nondimensionalization 
[6].) Previous experience with this flow [11] indicates the flow attains a statistically steady 
state more rapidly as the forcing wavenumber is increased; the use here of kf = 6 provides 
fairly brief initial transients. 

In the direct numerical simulation that provides the data set employed in the a priori 
analysis, the code of Shebalin [12] was used as configured for the earlier DNS found in 
Shebalin and Woodruff [7]. This code is pseudo-spectral with Fourier modes in the three 
spatial directions. A predictor-corrector time stepping algorithm is employed, with the 
viscous terms handled implicitly and the nonlinear terms handled explicitly. 

The large-eddy simulations were performed with the same code, as adapted for LES 
for the simulations discussed in Woodruff, Shebalin and Hussaini [11]. The sub-grid-scale 
stresses are incorporated as an explicit contribution to the time stepping; in this case the 
traditional Smagorinsky expression is used [13] 


T~i j — C5 A \jS mn S mn Sij , 


( 2 ) 
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where C s is the Smagorinsky constant, A 2 is a measure of the grid size and Sij is the rate- 
of-strain tensor. As has been customary with LES using the Smagorinsky model and Fourier 
modes, the grid size (relative to the length of the sides of the box) is taken as 2n/N, with 
N the number of modes in each coordinate direction. (All simulations considered here have 
the same number of modes in each coordinate direction.) 

Since the fast-Fourier transform algorithm in the code expects the number of modes to 
be a power of two, runs with intermediate numbers of modes are performed by zeroing out 
the “excess” modes at each time step. For example, 24 3 runs are performed with the code 
set for 32 3 and with modes with wavenumbers of 13 through 16 zeroed out. 

All DNS and LES runs reported here were initialized with a random initial velocity field 
with an exponential energy spectrum. Time series for the total kinetic energy, dissipation, 
etc. were used to establish that a statistically steady state had been reached. 

A priori Analysis. 

The goal of the a priori analysis is to determine what grid-size dependence gives the 
best fit between the sub-grid-scale turbulent stresses predicted by the Smagorinsky model 
and the true sub-grid-scale stresses when both are determined from a DNS data set. A 64 3 
simulation was run as described in the previous section until a statistically steady state was 
reached, then the instantaneous velocity field from this simulation was used in the analysis. 

The instantaneous velocity field was first filtered into super- and sub-grid components 
by means of a spectral cut-off filter (denoted (•)), which is to say that the velocity field was 
decomposed according to 

u i = ( u % ) + u 'i ( 3 ) 

with (u^ containing all Fourier modes with wavenumber less than the cut off k c and u\ = 
u . _ { u f ) containing all Fourier modes with wavenumber greater than k c . This decomposition 
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■rformed for all integer values of k c running from 32 (no filtering, with all modes 
led in the super-grid velocity u t ) down to 2 (nearly complete filtering, with only the 
k = 1 modes contained in the super-grid velocity and all the others contained in 

b-grid velocity u'.) 

each of the cut-off wavenumbers k c , the Smagorinsky representation of the sub-gnd- 
.tresses was computed, as well as the true stresses it is supposed to predict. The two 
s An = Sij and r 0 - = - <«$«}) that result are, according to the Smagorinsky 

(eqn. 2), supposed to be proportional, with proportionality factor C s A 2 . An empirical 
of the proportionality factor may be determined by reducing these second-rank tensor 
varying in x, y and z to single numbers and taking their ratio. 

ch of the components of these two tensors is a function of x, y and 2 ; in this work, 
inction is reduced to a single number by taking the L 2 norm. Such a norm seems 
priate for this problem, given the statistical homogeneity of the flow in the x and y 
ions. The flow is not statistically homogeneous in the 2 direction, but even there flow 
es do not cluster in particular regions of the flow domain and a norm which weights 
y all parts of the flow seems reasonable. A more specialized norm will probably be 
;ary in order to get useful results when such an analysis is applied to, say, turbulent 
lary-laycr flow or turbulent channel flow. 

summary, then, for each filter cut-off wavenumber the velocity field is decomposed 
upergrid and subgrid parts, the true sub-grid-scale stress and the Smagorinsky sub- 
tress are formed, the L 2 norms are taken of each of the six independent components 
. tensors, and the ratios are taken. In the results presented here, these ratios have 
further divided by A 2 , so that the result of the calculation is a grid-size dependent 

orinsky constant Cs- 


does exist. However, at this early stage of this research program, this measure of agreement 
provides a simple and surprisingly effective means of comparison. 

In addition to the 32 3 reference simulations, large-eddy simulations were carried out at 
resolutions of 24 3 , 16 3 , 14 3 and 12 3 . Time histories of results are shown in Figures 2-6; 
Table 1 contains parameters for the individual simulations and time averages of representative 
quantities over the statistically steady portion of the simulations’ time histories. The time 
histories of the total kinetic energy for the most extreme case, the 12 3 simulations, are 
compared against the 32 3 reference simulations in Figure 2. Time histories are shown for 
simulations with values of the Smagorinsky constant that bracket the 32 3 reference plots; the 
two time histories are shown to give a rough impression of the sensitivity of the results to 
the change in Cs at each resolution. Thus, at least as far as the kinetic energy is concerned, 
adjustment of the Smagorinsky constant permits successful simulations for resolutions as low 
as 12 3 . The extent to which other quantities are successfully simulated is examined in detail 
below. 

Interpolated values for Cs were plotted in Figure 1 for comparison with the Cs values 
determined from the a priori analysis; it is seen that the LES experiments are consistent with 
the a priori results. The determination of similar variations in the Smagorinsky constant 
by these two independent means is strong evidence for the utility of alternative grid-size 
dependences in LES. Specifically, these results indicate that not only is the enhanced grid- 
size dependence necessary for low-resolutioned LES (as indicated by the a priori analysis), 
but it is sufficient (at least as far as the kinetic energy is concerned.) 

Time histories of some other first and second order statistical quantities are shown in 
Figures 3-6. In each Figure, one simulation at each resolution is shown. (In all cases, the 
simulation is the one with the higher kinetic energy.) The statistical quantities are computed 
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using averages in the x - y plane. Thus, u is a function of z; Figure 3 shows the L 2 norms 
of these functions. The fluctuating velocities u' — u- u and w' = w — w are represented by 
their L 2 norms in Figures 4 and 5. (Note these fluctuating velocities are different from the 
primed, fluctuating, velocities defined in the previous section.) The final plot, Figure 6, are 
of the L 2 norms of the correlation coefficient, defined by 



Considering first the mean and fluctuating velocities of Figures 3-5, it is clear that 
almost all the simulations’ predictions for these quantities cluster together quite well, once 
the initial transients have evolved into the statistical steady state. The exception is the 
12 3 simulation, where it is clear that the limits of resolution have been exceeded. It is 
interesting to note that the 12 3 simulation has achieved the correct total kinetic energy 
by significantly overestimating the plane-averaged mean velocity and underestimating the 
fluctuating velocities, indicating the low resolution is suppressing the fluctuations. 

The time series of the L 2 norm of the correlation coefficient shown in Figure 6 provide a 
more rigorous measure of the faithfulness of the simulations to the true turbulence physics, 
and their examination reveals weaknesses in the lower-resolutioned simulations. The C xy 
and C yz correlations (not shown) are reproduced fairly well by all simulations (even the 12 3 
simulation); it is, however, the C xz correlation, describing the Reynolds stress component 
that effects energy transfer from the mean to the fluctuating motion, where real breakdowns 
in the low-resolution simulations may be seen. Figure 6 shows that only the 24 3 simulation 
reproduces the steady-state value of C xz correctly, with the 16 3 and 14 3 simulations giving 
values around half the true value. Curiously, and probably fortuitously, the 12 3 simulation 
does better here than the 14 3 and 16 3 simulations. It is worth noting as well (see Table 1) 
that the correlation coefficients are essentially unaffected by the changes in the Smagorinsky 
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constant at a given resolution. 


Summary. 

These results indicate that the low-resolution simulations not only correctly reproduce the 
total kinetic energy, but, except for the lowest-resolution, 12 3 simulation, correctly reproduce 
most of the other statistical quantities examined as well. The major failure of the low- 
resolution simulations was in the determination of the correlation coefficient C xz , where only 
the 24 3 simulation came close to the correct value. 

In addition to this failure of the low-resolution LES results in the prediction of the 
C xz correlation, the results (as presented in Table 1) also show virtually no sensitivity of 
this correlation coefficient to changes in the value of the Smagorinsky constant. It may be 
concluded, then, that not only do the simulations performed with the Smagorinsky constant 
tuned to give the correct kinetic energy fail to give the right C xz , but that no simple tuning 
of the Smagorinsky constant will give good results for the correlation coefficient. A more 
drastic adjustment of the model is necessary, and is being developed. 

The breakdown of the under-resolved simulations at 12 3 and the failure of the under- 
resolved simulations to produce the correct correlation coefficient C xz are likely related. 
Reducing the number of modes and the resolved scales inhibits the mechanism for the transfer 
of the energy from the mean to the fluctuating motion. At the 12 3 resolution, this inhibition 
has reached the point where insufficient energy flows into the fluctuating motion and so the 
mean energy is too high and the fluctuating energy is too low. Since, in this flow, as in other 
parallel flows, the energy transfer between the mean and fluctuating motion is mediated 
by the C xz correlation, the inhibited energy transfer is reflected in incorrect values for that 
correlation. 

The a priori results and the LES experiments described here indicate that, when used 
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with care, low-resolution simulations with an alternative grid-size dependence can give sat- 
isfactory results. Significant advances in the application of LES to complex turbulent flows 
of practical interest may thus be possible, provided reliable, general-purpose models are de- 
veloped incorporating the concept of alternative grid-size dependences. Research into such 
models is currently under way. 
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Table 1: Simulation parameters and temporal averages of results. 


No. Modes 


K 

INI 

Ill'll 

lb'll 

\w' 

c zx 

32 

0.006 

2.33 

0.93 

1.24 

0.89 

1.21 

0.41 

32 

0.008 

2.08 

0.74 

1.20 

0.82 

1.22 

0.44 

24 

0.0045 

2.48 

0.88 

1.36 

0.96 

1.19 

0.40 

24 

0.0056 

2.33 

0.85 

1.32 

0.90 

1.17 

0.41 

16 

0.008 

2.49 

0.96 

1.28 

0.88 

1.28 

0.25 

16 

0.012 

2.26 

1.04 

1.21 

0.71 

1.20 

0.27 

14 

0.0077 

2.38 

0.85 

1.26 

0.82 

1.33 

0.27 

14 

0.010 

2.18 

0.82 

1.15 

0.69 

1.37 

0.28 

12 

0.0087 

2.44 

1.75 

0.97 

0.81 

0.45 

0.32 

12 

0.0096 

2.25 

1.69 

0.92 

0.78 

0.42 

0.34 
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Figure 1: Smagorinsky constant Cs as function of grid-size, as determined by a priori analysis 
(solid line) and LES experiments (dots). 

Figure 2: Kinetic energy versus time for 32 3 simulation with Cs = 0.006 (solid line), 32 3 simu- 
lation with Cs = 0.008 (long-dashes), 12 3 simulation with Cs = 0.0087 (dashed line), 
12 3 simulation with Cs = 0.0096 (dotted line). 

Figure 3: ||u|| versus time for 32 3 simulation with Cs — 0.006 (solid line), 24 3 simulation with 
C s = 0.0045 (long dashes), 16 3 simulation with Cs = 0.008 (dashed line), 14 3 simu- 
lation with Cs = 0.010 (dash-dotted line), 12 3 simulation with C s = 0.0087 (dotted 
line). 

Figure 4: ||u , || versus time for 32 3 simulation with Cs — 0.006 (solid line), 24 3 simulation with 
C s = 0.0045 (long dashes), 16 3 simulation with Cs — 0.008 (dashed line), 14 3 simu- 
lation with Cs = 0.010 (dash-dotted line), 12 3 simulation with Cs = 0.0087 (dotted 
line). 

Figure 5: ||u/|| versus time for 32 3 simulation with Cs — 0.006 (solid line), 24 3 simulation with 
C s = 0.0045 (long dashes), 16 3 simulation with Cs = 0.008 (dashed line), 14 3 simu- 
lation with Cs = 0.010 (dash-dotted line), 12 3 simulation with Cs = 0.0087 (dotted 
line). 

Figure 6: ||C I2 || versus time for 32 3 simulation with Cs = 0.006 (solid line), 24 3 simulation 
with Cs = 0.0045 (long dashes), 16 3 simulation with Cs = 0.008 (dashed line), 14 3 
simulation with Cs — 0.010 (dash-dotted line), 12 3 simulation with Cs = 0.0087 
(dotted line). 
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Figure 2: Kinetic energy versus time for 32 3 simulation with C s = 0.006 (solid line), 32 3 
simulation with C s = 0.008 (long-dashes), 12 3 simulation with C s = 0.0087 (dashed line), 
12 3 simulation with C s = 0.0096 (dotted line). 
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Figure 3: ||u|| versus time for 32 3 simulation with C s = 0.006 (solid line), 24 3 simulation with 
C s = 0.0045 (long dashes), 16 3 simulation with C s = 0.008 (dashed line), 14 3 simulation 
with C s = 0.010 (dash-dotted line), 12 3 simulation with C s = 0.0087 (dotted line). 
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Figure 4: ||ti'|| versus time for 32 3 simulation with C s = 0.006 (solid line), 24 3 simulation 
with C s = 0.0045 (long dashes), 16 3 simulation with C s = 0.008 (dashed line), 14 3 simulation 
with C s = 0.010 (dash-dotted line), 12 3 simulation with C s = 0.0087 (dotted line). 
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Figure 5: ||tt/|| versus time for 32 3 simulation with C s = 0.006 (solid line), 24 3 simulation 
with C s = 0.0045 (long dashes), 16 3 simulation with C s = 0.008 (dashed line), 14 3 simulation 
with Cs = 0.010 (dash-dotted line), 12 3 simulation with Cs = 0.0087 (dotted line). 



Figure 6: ||C XZ || versus time for 32 3 simulation with Cs = 0.006 (solid line), 24 3 simulation 
with Cs = 0.0045 (long dashes), 16 3 simulation with Cs = 0.008 (dashed line), 14 3 simulation 
with Cs = 0.010 (dash-dotted line), 12 3 simulation with Cs = 0.0087 (dotted line). 
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